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, We study the behavior of the energy fluctuations in the stationary state of a uniformly heated 

granular gas. The equation for the one-time two-particle correlation function is derived and the 
hydrodynamic eigenvalues are identified. Explicit predictions are subsequently worked out for energy 
fluctuations. The results explain Monte Carlo numerical data reported in previous studies (P. Visco 
f-^ . et al, European Physical Journal B 51, 377 (2006)). 
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I. INTRODUCTION 



Recent years have witnessed ongoing interest for the microscopic and macroscopic properties of granular media. 
In such systems, a simple ingredient -energy dissipation resulting from collisions- has far reaching consequences 
[l|, with rich phenomenology: non-Gaussian velocity distributions @, Q, non-equipartition of energy 0, [ESQ, ° r 
spontaneous symmetry breaking 0, [13] to name but a few. Theoretically, one of the tools used to understand 
this body of phenomena is kinetic theory, which is extended naturally to these systems by introducing an "inelastic 
collision rule" in which the energy is not conserved. Although most of the work carried out until now has focused 
on the one-particle properties, and on the study of the corresponding Boltzmann equation, it has been shown that 

S correlations are also important and, as a matter of fact, necessary to understand the behavior of the system when 
vortices or cluster are developed [ll|, [l2j , or even in simpler situations where spatial homogeneity is enforced (l3l . [lij . 
, As a consequence of the dissipation in collisions, the total energy of an isolated granular system decays monotonically 
in time. Under certain conditions, the system reaches a homogeneous cooling state in which the time dependence of the 
one-particle distribution function is entirely embodied in the kinetic (so-called granular) temperature, which evolves 
with time as t~ 2 It is experimentally difficult to probe such a regime (see however [III 113] )• Nevertheless, it 

is possible to maintain a granular system in the fast-flow regime by injecting energy in such a way that a stationary 
y—i , state is reached. In these states, the energy injected by the thermostat is compensated by the energy dissipated in 
' collisions. Several mechanisms can be introduced in order to get a stationary state. If, for example, energy is injected 
by a moving boundary such as a vibrating piston, the system reaches an inhomogeneous stationary state In this 
work, we will focus on a granular gas which is heated uniformly, coupling the velocity of each particle to a white noise, 
this is the so-called "stochastic thermostat" 0, H3, HU, [22|, [H, |2l, [H, H3] - For this kind of forcing, which is 
relevant for some two-dimensional experimental configurations with a rough vibrating piston [28| , the system reaches 
a homogeneous stationary state after a transient regime. The advantage of such a driving mechanism is that it lends 
itself to theoretical progress. In this context, the single particle distribution function has been characterised 0] and 
long range correlations predicted from a Jiydrodynamic treatment [2l|. More recently, the fluctuations of the total 
energy have been analyzed [29[ (see also [1q| for a related numerical study in an inhomogeneous system). In Ref. [29| . 
the second moment of the total energy fluctuations was evaluated by neglecting the correlations, which, by and large, 
could not explain the simulation results. The objective of the present work is to clarify and quantify the influence of 
the inelasticity induced correlations on the total energy fluctuations. The methods used bear some similarities with 
those reported in [3l| , where it was shown that for the unforced system, the contribution coming from the correlations 
is of the same order as that coming from the one-particle distribution function itself. 

The paper is organized as follows. In section [TlJ the equation for the two-particle distribution function is derived, 
taking due account of the thermostat while in section IHT1 the results are particularized to the homogeneous stationary 
state, that will play the role of our reference state in subsequent analysis. There we also summarize the main results 
already known and pertaining to the one-particle distribution function. In section HVi we analyze the hydrodynamic 
equations for a homogeneous linear perturbation of the reference state, and obtain the corresponding modes and 
eigenvalues, which are finally used in section [V] to meet our objective and obtain an explicit expression for the 
variance of the total energy. 
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II. HEATED GRANULAR GAS : TWO BODY KINETIC DESCRIPTION 



We consider a gas of N hard disks (dimension d = 2) or spheres (d = 3) of mass m and diameter a that collide 
inelastically with a coefficient of normal restitution a The system is heated uniformly by adding a random 
component to the velocity of each particle at equal times [U [12] ■ The driving is implemented in such a way that 
the time between random kicks is small compared to the mean free time. Then, between collisions, the velocities of 
the particles undergo a large number of kicks due to the thermostat. In addition, we will assume that the "jump 
moments" of the velocities of the particles verify 

i,j = l,...,N and (3,^ = 1, . . . ,d 

where we have introduced Avi^ = Uj p(t + At) — Vi p(t), Vi being the (3 component of the velocity of particle 
i at time t. We have also introduced the strength of the noise, £qi an d (•••)) which denotes average over different 
realizations of the noise. The non-diagonal terms (corresponding to i ^ j and (3 = 7) are necessary in order to 
conserve the total momentum. 

In the dilute limit, assuming molecular chaos, i.e. that no correlations exist between colliding particles, and that the 
sizes of the jumps due to the thermostat are small compared to the velocity scale on which the distribution varies, the 
equation for the single particle distribution function in our system is the Boltzmann-Fokker-Planck equation [3|, |22| | 

— f(T., t\ 4- rSV(Ti)f(T.,.f) = j\f\n + ^ f —\ 



-f( Xl ,t)+L^\x 1 )f(x 1 ,t) = J[f\f} + ^l — \ f( Xl ,t), (2) 



where short-hand for position- momenta coordinates {r^Vi} and 

_d_ 

dri 



i (0) (o :i )=v 1 .A. (3 ) 



The inelastic collision operator J[/|/] reads 

J[f\f] = <r d - 1 J dv 2 T (v 1 ,v 2 )/ 1 (r,v 1 ,i)/ 1 (r,v 2 ,t), (4) 

where 

f„(vi,v 2 ) - J d&0(& ■ g)(«r • g) [a-% 1 - 1] , (5) 

with g = vi — v 2 the relative velocity, the Heaviside step function, <x a unit vector joining the centers of the particles 
at contact and b^. 1 an operator replacing the velocities vi and v 2 appearing on its right by the precollisional values 

v* = fc-'v! = vi - ^^(g- (6) 
1 + Oi 

V 2 = b » lv 2 = V 2 + (g ■ <T)<T. (7) 

The term ^ (^§7) f( x i>t) 1S a diffusive Fokker-Plank term, and is a signature of the external noise. 

As we shall study fluctuations, it is convenient to introduce the two-particle distribution function, f2(x\,X2,t). 
The quantity fi{xi,X2,t)dxidx2 is defined as the number of pairs of particles in which one lies inside the differential 
volume dx\ centred in x\ and likewise, with dx2,X2 for the second particle. This definition is easily generalized to 
higher n-particle distribution functions, f n (xi, ■ ■ ■ , x n )- The evolution equation for / 2 is [3l|, [32j 

|-+£ (0) (£l) + £ (0) (z 2 ) f 2 (xi,X2,t) = 6(T 12 )<T d - 1 f (v U V2)f2(xi,X2 > t) 

at 

+<r d ~ 1 [ dx 3 [5(r 13 )T ( Vl ,v 3 ) + <5(r 23 )r (v 2 ,v 3 )] f 3 (xi,x 2 ,x 3 ,t) + F TH , (8) 

where we have introduced Fth that accounts for the external driving. The evolution equation ([5]) contains essentially 
three parts: the free streaming in the left-hand side, the two terms in the right hand side corresponding to collisions, 
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and the last term, Fth due to the thermostat. The collisional contribution is split in one part corresponding to 
collisions of particles with velocities Vi and v 2 , and the other which involves collisions of particles with velocities Vj 
or V2 with a third particle with arbitrary velocity, V3. The collisional contribution is identical to the one that appears 
in the absence of forcing (3l| . We concentrate now on the new term, Fth- Assuming that the sizes of the jumps 
due to the thermostat are small compared to the scale in which the distribution f 2 varies, we can expand Fth i n the 
spirit of the Fokker-Planck description [33| 



Fth 



1 d 2 

\ E E % 

0,7=1 «,j=l 



9 9 ff A 



2 to 



dv hl3 dv jn 
2 d 



d 2 



(9) 



where we have taken into account equation l[T]). and we have explicitly assumed that the jump moments By do not 
depend on the magnitude of the velocities of the particles. 

Let us introduce the two-particle and three-particle correlation functions through the usual cluster expansion 



(10) 



and 



h{xux 2 ,t) = f 1 (x 1 ,t)f 1 (x 2 ,t) + g 2 (xi, x 2 ,t), 

f 3 (x 1 ,x 2 ,x 3 ,t) = fi(x 1 ,t)f 1 (x 2 ,t)f 1 (x 3 ,t) + g 2 (x 1 ,x 2 ,t)f 1 (x 3 ,t) 
+92 (xi,x 3 , t)ji (x 2 ,t)+g 2 (x 2 , x 3 , t)fi (xi , t) + g 3 {x x , x 2 , x 3 , t) . 

The equation for the correlation function g 2 (xi, x 2 , t) can be obtained following the same lines as in references [3ll. [34 
Neglecting the three-body correlations, g 3 , in Eq. ©, we obtain 



(11) 



^-+L^\x 1 )+L^{x 2 ) 



g 2 (xi,x 2 ,t) = 8(r 12 )(j d 1 T (v 1 ,v 2 )f 1 (x 1 ,t)fi(x 2 ,t) 



+ \K{x u t) + K(x 2 ,t)} g 2 { Xl ,x 2 ,t) - • 7^/1 (zi, t)h(x 2 , t), 



where we have introduced the linear operator K(xi, t) defined as 
K( Xi ,t) = a d ~ 1 



dx3<5(r l3 )f (v 4 ,V3)(l +V i3 )fi(x 3 ,t) + f f A 

2 \dvi 



(12) 



(13) 



and where the permutation operator V a b interchanges the labels of particles a and b in the quantities on which it acts. 
As will become clear below, the 1/N term in equation (fT2|) is crucial for the calculation of the energy fluctuations. 



III. THE STATIONARY STATE 

It has been shown numerically that, after a transient time, the system reaches a homogeneous stationary state [2l| 
in which the energy input from the thermostat is compensated by the energy lost in collisions. In this section we 
will particularize the equations of the previous section to this state, summarizing the results that are already known 
about the one-particle distribution function and that are required for our theoretical analysis. 

The Boltzmann- Fokker-Planck equation ([2]) for the distribution function, /#(vi), in the stationary homogeneous 
state is 



^-(J-XfH^l) + J[fH\fH}=0. 



2 V<9viy 

It is convenient to introduce the scaled distribution function \h 



(14) 



Wv) = %W, (15) 



where nn is the homogeneous density, vh = (^ L ) 1 ^ 2 is the thermal velocity defined from the granular temperature 
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and c = ^p- is the rescaled velocity. The distribution function has been studied in reference [3| , where an approximate 
expression for xh(c) was derived to second order in Sonine polynomials (35j 

X"( C ) = ( l + ^i a )S 2 d /2-l(c 2 )) , (17) 



with 

S^/a-iCc 3 ) = 5c 4 - i(d + 2)c 2 + id(d + 2), 
and 02(0;) a coefficient related to the kurtosis of the function xh(c) 

d {c A ) H 



2 (c 2 ) 2 H 



= 1 + 02(a). 

An approximate expression for a 2 reads (see [23l l36j | for a discussion on various possible approximations) 

_ 16(1 - a){\ - 2a 2 ) 

° 2 ~~ 73 + 56d - 24ad - 105a + 30(1 - a)a 2 ' 

The expression for the temperature in the first Sonine approximation is 



(1 - a 2 )n d n H a 



d-l 



(l + 0(a 2 )), 



(18) 



(19) 



(20) 



(21) 



where fi<j = 2n d ^ 2 /T(d/2) is the d-dimensional solid angle. 

We now turn to the equation for the correlation function, g2,H(%i, #2)- It is convenient to introduce the rescaled 
correlation function gn through 



g2,H{xi,x 2 ) 



e d vjf 



gH(h2,Cl,C 2 ), 



(22) 



where I — (nno d v ) 
function gn reads 



_1 is proportional to the mean free path and 1 = r/£. In these units, the equation for the reduced 



d 



A(ci) + A(c 2 ) - C12 • — — 
CI12 



3_f/(ll2,Ci,C 2 ) 

S(h2)T (c 1 ,c 2 )x H (ci)XH(c2) + £0 



- 2 n H e d d 
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N dci dc 2 

where we have introduced the linearized Boltzmann-Fokker-Planck operator A(c) 



Xh(ci)xh(c 2 ) 



V(C, )MC i = / dc 3 r (c,;, Cs)(l + P t 3)XH(c3)h( Cl ) + ^ 



d_ 

dci 



with rescaled noise amplitude 



Co 



v 3 H 



(23) 



(24) 



(25) 



As can be seen in equation (|23|) . the correlation function, gn, is determined by the properties of the linearized 
Boltzmann-Fokker-Planck operator, A, and by the one particle distribution function xh- It is consequently important 
to study the spectral properties of A, in particular the upper (hydrodynamic) part of the spectrum, in order to 
understand the fluctuations of global quantities. In the case of a granular gas in the homogeneous cooling state (2Z[, and 
for a system under ballistic annihilation dynamics [381 ] , it has been shown that it is possible to find the hydrodynamic 
eigenvalues and eigenfunctions of the linearized Boltzmann-Fokker-Planck operator. Once these quantities are known, 
it becomes possible to evaluate the fluctuations of the relevant global quantities in the so-called "hydrodynamic 
approximation" [3ll . [34| . In the remainder, we will see that we can evaluate the fluctuations of the total energy in 
an equivalent approximation, but without the knowledge of the eigenfunction associated to the energy. The only 
information needed is the form of the linearized hydrodynamic equations and, in particular, the eigenvalues. 
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IV. HYDRODYNAMIC EQUATIONS 
A. Evolution of homogeneous perturbations 

In this section we focus on the linearized hydrodynamic equations around the homogeneous stationary state. The 
objective is to consider the linearized equations around a homogeneous perturbation in order to extract information 
about the linear behavior of a small perturbation of the total energy. 

The complete non-linear hydrodynamic equations for the granular system heated by the stochastic thermostat are 

[mm 

^-n = -V-(nu), (26) 



-|u = -u- Vu- — VjPij, (27) 
at mn 



jT = -u ■ VT - J^(V • q + PtjVjUi) - CT + (28) 

where i^y is the pressure tensor, q is the heat flux and £ is the cooling rate, which is also a functional of the distribution 
function 

( - / dTl / " V2|T1 " V2|3/(r ' T " ' )/(r ' Y2 ' (29) 

Considering a homogeneous state, the previous equations reduce to 

|» = 0, |u = 0, §-T = -(T + m e . (30) 

In the long time limit, the system is expected to approach a steady state with a constant temperature given by the 
equation 

(H(f H )T H =mZl (31) 



Substituting the explicit form of the one particle distribution function (|15[) in the equation above, we obtain the 
temperature given in equation (f2"Tj) . 

Let us consider now a homogeneous state close to this homogeneous stationary state. We can write the hydrodynamic 
fields as n(t) = njj + fin, u(t) = Su and T(t) = Th + ST. We also define the dimensionless hydrodynamic fields 

Sp(r) = —, Sw(t) = —, 59(t) = ^, (32) 
n H v H 1 H 

where we have introduced the dimensionless time scale r, proportional to the number of collisions per particle, defined 

as 

V^ = ft. (33) 

Assuming that the deviations are small, and taking into account equations (|30p - (|31|) . we can write the linearized 
evolution equations for the dimensionless hydrodynamic fields in this new time scale 

^s P = o, J^w = o, ^se = -CoS P -^se, (34) 



VH 

stationary state. Its expression in the first Sonine approximation is [3J 



where Co = is a dimensionless coefficient that is a functional of the one-particle distribution function in the 



(16 + 3^)7^(1 -a 2 ) 

Co = sTl^ • (35) 

To obtain the equation for SO we have assumed that the perturbed distribution function scales as 



6 



where v(t) — , and xh is the same scaled distribution function as for the reference stationary state. This 

assumption has already been used and tested numerically in (2l| . Then, the cooling rate C for the state under scrutiny 
is proportional to T 1 / 2 (i) and we obtain the equation for the linearized energy written above. Equations (|34j) indicate 
that a perturbation in the total number of particles or total momentum does not decay, as a consequence of the 
fact that these variables are conserved, but a perturbation in the total energy will decay (exponentially in r) to the 
stationary value, as expected. Moreover, as the equation for the temperature can be rewritten in the following form 

L(ls p + 86\=-l<j>(\5p + 69\, (37) 



1/2 



dr \3 ) 2^ u V 3 

we can identify the hydrodynamic eigenvalues A = and 7 = — f Co, A being [d + l)-fold degenerate. For the sake of 
clarity, it proves convenient to relabel these eigenvalues as 

A 1= 0, A 2 = 0, A 3 --^Co, (38) 

where A2 is G?-fold degenerate. The associated hydrodynamic modes, {yp}pt?i are 

2 

2/i = Sp, y 2 = £w, 2/3 = -^Sp + 50. (39) 

The + sign in the last equation stems from the fact that an increased density leads to enhanced dissipation, and 
hence, a smaller temperature. 



B. Enforcing consistency with the linearized Boltzmann-Fokker-Planck equation description 

We now turn our attention to the problem of finding the linearized hydrodynamic equations for a homogeneous 
perturbation, directly from the Boltzmann-Fokker-Planck equation. Enforcing consistency with the macroscopic 
considerations of section IIV Al above, we will infer useful properties on the hydrodynamic part of the spectrum of 
A(c). We first introduce the scaled deviation of the distribution function 

5x(c,r) = ^[/(v,t)-/ ff (v)]. (40) 
n H 

The evolution of the scaled distribution is governed by 

JUx(cr)=A(c)*x(cT), (41) 

where the operator A(c) is the linearized Boltzmann-Fokker-Planck operator defined in (|24[) . Let us also introduce 
the scalar product 



(/(c)b(c)) = J de X x L (c)r(c)g(c), (42) 

where /* is the complex conjugate of /. Interestingly, the hydrodynamic modes introduced in (I39|) can then be written 
as 

VP = &X), = 1,2,3, (43) 

where 

6(c) = Xh(c), 6(c) - c XH (c), 6(c) = (7 - ^) Xh(c). (44) 

Taking the scalar product of the linearized Boltzmann-Fokker-Planck equation (|4"Tj) with the functions we obtain 
the linear equations (j3~4"|) (in the hydrodynamic time scale, that is, if we wait long enough so that fast modes have 
vanished) only if the spectrum of A admits the three eigenvalues written in ([38]) . and the associated "hydrodynamic" 
eigenf unctions, {^/slsii) °bey the orthogonality condition 



(£81 l£ft, ) = <W 2 , 01) A = 1,2,3. 



(45) 
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In Appendix [5] it is shown that the null eigenvalue is (d + l)-fold degenerate, and the corresponding eigenfunctions, 
£1 and £2 > are worked out. Moreover, as a consequence of particle and total momentum conservation in a collision, 
£1 and £2 are the corresponding left eigenfunctions. We were not able to demonstrate that A3 is an eigenvalue of A, 
but we have shown explicitly that 

(&|&> = 0, for 0=1,2. (46) 

In the following, we will assume that A actually admits this third eigenvalue, with an unknown eigenfunction £ 3 . With 
the help of this assumption, we will see in the next section that it is possible to define a projector in the hydrodynamic 
subspace, which opens the way for evaluating the variance of the total energy fluctuations. 



V. ENERGY FLUCTUATIONS 



In this section, we study the fluctuations of the global energy for a system in the stationary state. As we are 
interested in global quantities, it is convenient to define a global correlation function </>jj 

0ff(ci,c 2 ) = J dri23ij(ri2,ci,c 2 ). (47) 

The energy fluctuations can be written as a functional of this correlation function and the one-particle distribution 
function \h, as 



{{5Ef) H = ^Nv% 



J dcc 4 XH(c) + j ' dci J dc 2 cl^ H (ci,c 2 ) 



(48) 



In order to evaluate the integral over the correlation function we start from (|23[) . integrating over the position 
variable. Assuming periodic boundary conditions, the spacial gradient terms disappear and we have the following 
equation for 

[A(ci)+A(c2)]Mci,C2)=r(ci,C2), ( 49 ) 

where 

~ d d 

r(ci,c 2 ) = -T (ci,c 2 )Xh(ci)xh(c2) +£o7^r~ 1 -q^Xh{ci)xh{c 2 ). (50) 

The solubility condition for equation (|49p is that T does not have components in the subspace associated to the null 
eigenvalue. In our case, this subspace is generated by {£i)£a}. Due to the conservation of the number of particles 
and total momentum in a collision, and to the symmetry of the second term of T, we have 

{ii{ci)^(c 2 )\T( Cl , c 2 )) = (f 2 (c 1 Ki(c 2 )|r(c 1 ,c 2 )> = 0, (51) 
(a(ci)6(c 2 )|r( Cl ,C2)) = (^A^)kj{^)\T(c u c 2 )) = 0, (52) 

for i 7^ j. The case i = j is analyzed in the Appendix [Bl where it is shown that (£ 2 ,i(ci)£ 2) j(c 2 )|r(ci, c 2 )) = 0. In 
order to prove this property, the presence of the second term of the right hand side of Eq. ([50)) is essential. Hence, 
the solubility condition holds and the problem of finding (f>H with equation (|49)l is well defined. 
Let us also define a projector P\ 2 in the hydrodynamic subspace as 

3 3 

Pi 2 %i,c 2 ) = ^2 ^2(^l( C l)^2( C 2)\Hci,C2))^} 1 {c 1 )C f 3 2 {c 2 ) 1 (53) 

01=1 03=1 

where {£/3}j| =1 are the right hydrodynamic eigenfunctions of the linearized Boltzmann-Fokker-Planck operator, 
{£/3}fl=i tne orthogonal set introduced in the previous section, equation (|44]) . and we have generalized the scalar 
product by 

(/(ci,c 2 )|g(ci,c 2 )) = J da J dc2X ff 1 (ci)x ff 1 (c 2 )/*(ci,c 2 )g(ci,C2). (54) 

Note that Pi 2 is a projector even if the function £3 is not the true left eigenfunction of A(c) (remember that £1 and 
£ 2 are the actual left eigenfunctions associated to the null eigenvalue). The fact that the set of functions and 
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{f/3}jg = i fulfil the orthogonality condition (|45|) is enough to guarantee that Pf 2 = Pi 2 . Using this projector, we define 
the "hydrodynamic part" of 4>h to be the function 

3 3 

0h (ci,c 2 ) = Pi 2 0h(ci,c 2 ) = ^2 a fJi^f3 1 (ci)^ 2 (c 2 ). (55) 

01 = 1 02=1 

The coefficients ap x p 2 are the quantities we need to evaluate. As they are essentially the first moments of the correlation 
function 4>h, they are directly related to the integral we have to calculate in (|4"5| . It is tempting to treat £3 as if it 
was an actual left eigenfunction of the linearized Boltzmann-Fokker-Planck operator, and we will in the following use 
the approximation 

P 12 A(c,) = P 12 A( Ci )P 12 , (56) 

which allows to find a closed equation for . This approximation has already been invoked in other systems such 
as the freely evolving granular gas [3l[, or t ne probabilistic ballistic annihilation model [HJ|. With the information 
available on the linearized Boltzmann-Fokker-Planck operator, it seems the best that can be done technically. Let us 
also remark that the approximation is exact in the elastic limit. Then, applying the projector Pi 2 to the equation 
(|49|) and taking into account the approximation (|56p . we obtain the following expressions for the coefficients 

(6l( C l)Cfe( C 2)|Ji)(ci,C2)x/j-(ci)xg(c2)) , . 

Wife = r — — t , (57) 

where it has been assumed that \p t + Xp 2 7^ 0. The coefficients associated to the vanishing eigenvalue cannot be 
calculated by equation (|49p . but are fixed by the boundary conditions. In Appendix [Cl the coefficients ap r p 2 are 

evaluated. The expression for <j>^ is finally given by 

4$ = oil€i(ci)6(c 2 ) + 013 [6(ci)6(c 2 ) + 6(ci)a(c 2 )] + a 3 36(ci)6(c2), (58) 
where an = —1, ai3 = — i and 033 can be obtained as a functional of the one-particle distribution function 

(&(ci)6(c 2 )|IL(ci,c 2 )xff(ci)x H (c 2 )) 

a 33 57 • (5yj 

•Jt.0 

An approximate expression is derived in appendix [C] and reads 

_ -15 + 7d+ Ud 2 - 3(-9 + d(9 + 2d))a + 30(1 + d)a 2 - 6(9 + d)a 3 ) 
° 33 _ 9<f(— 19 + 2d(-7 + 3a) + 3a{9 + 2(-l + a)a)) ' ( ' 

Taking into account and the variance of the energy fluctuations can finally be written as 

E= (Pp = ( a 2 + l) 1 + da 3 3-d^- (61) 

In reference [13 > the value of cr^ has been measured by means of the Direct Monte Carlo simulation method (DSMC). In 
Fig. [I] we compare our theoretical prediction (solid line given by (I61| ) with the DMSC simulations results (symbols). 
The agreement is satisfactory for the whole range of inelasticities, at variance with the theoretical attempt put 
forward in [29j . which neglected velocity correlations. In particular, we note the non-trivial result in the elastic limit 
lim Q ^i- <7p(q) = d/3 (i.e. 2/3 on Fig. [TJ which is well obeyed), while in the free cooling regime, this quantity 
vanishes [3l|]. We emphasize that the elastic limit is singular: the behaviour for elastic systems with a = 1 is not 
approached by taking the quasi-elastic limit a — > l - (we note that the divergence of the different moments of the 
velocity distribution as a — > 1 is nevertheless indicative of the absence of a stationary state when a = 1). Such a 
singularity has already been reported in 1 dimension 39], but, to our knowledge, not for two dimensional granular 
systems. It is also interesting to note that the singular nature of the quasi-elastic limit appears at 2 body level through 
the energy fluctuations, while, as far as rescaled distribution functions are considered, the 1 body level of description 
is regular, with a well behaved velocity distribution approaching a Gaussian form [22j |. 
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Figure 1: Scaled second moment of the energy fluctuations cr% as a function of the restitution coefficient a. The solid line is 
the theoretical prediction and symbols are the two dimensional Monte Carlo simulation results of Ref. [29j ] . 

VI. CONCLUSIONS 

The problem of the fluctuations of the total energy of a granular (inelastic) gas maintained in a non-equilibrium 
stationary state by a random acceleration has been addressed. A numerical study of this quantity had been performed 
by means of Monte Carlo simulations and an argument assuming uncorrelated non-Gaussian individual distribution 
function had been proposed in [2{|, without success. The main goal of this work was therefore to take due account of 
velocity correlations in order to study these fluctuations. 

To this end, the standard description at the single particle level is not sufficient, and the two-particle correlation 
function is needed. We have derived the evolution equation for such an object, and particularized the analysis to the 
homogeneous stationary state that is reached by the system in the long time limit. Our work shows that this equation 
is not a straightforward generalization of its counterpart arising in the context of the undriven granular gas (i.e. by 
only changing the linearized Boltzmann-Fokker-Planck operator into its driven form). A non trivial non diagonal 
term appears in the Fokker-Planck equation for the two particle distribution function (contribution d Vl ■ d V2 in Eq. 
(|9"|)), as a consequence of the coupling between velocities due to momentum conservation. 

We have seen that, for our purposes, the exact knowledge of the hydrodynamic eigenfunctions is not needed. The 
important point is that we can construct a set of functions {£/3}j| =1 , which are linear combinations of 1, c and c 2 , 
that are orthogonal to the right eigenfunctions {^/3}^ =1 of the linearized Boltzmann-Fokker-Planck operator A. This 
orthogonality property holds for the "real" left eigenfunctions, that in our case correspond to the null eigenvalue (i.e. 
density and velocity fields associated to conserved quantities). The function £3 is not a left eigenfunction of A but it 
can be proved to be orthogonal to £1 and £2- In a subsequent step, the linear hydrodynamic equations around the 
reference state are derived and from that knowledge, the hydrodynamic eigenvalues are identified and the variance of 
energy fluctuations subsequently derived. 

Finally, our prediction has been successfully tested against the numerical results obtained by the Direct Monte 
Carlo simulation method for all the range of values of the coefficient of normal restitution a. This provides strong 
support for the theory developed here and assesses in retrospect the validity of our assumptions. 
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Appendix A: EIGENVALUE PROBLEM FOR A 

We consider here the eigenvalue problem for the homogeneous linear Boltzmann-Fokker-Planck operator A, defined 
in (22) 

A(c)&(c) = A^c). (Al) 

We are interested in the eigenf unctions and eigenvalues associated with linear hydrodynamics and, to perform the 
analysis, similar techniques as in [3l], [3?], [38| will be required. 
Consider first the function 

Mc) =Xh(c). (A2) 
When the linearized operator A acts on xh, we have 

A(cx) X h(ci) = J dc 2 f (c 2 ,c 3 )(l+V 12 )xH(c2)XH(ci) + ^ (J^\ xh(ci). (A3) 
Taking into account the equation for xh, Ec. (fl4f . we obtain the following relation 

A( Cl )Ma) = —(■£-) Xh{ci). (A4) 



2 V^ci 

Now, let us considerer the function 

Mc) = ~ik XH{c) - (A5) 

Taking dcrivate in the equation obeyed by xh(c — w) with respect to w, and subsequently evaluating the result for 
w = 0, we obtain 

A( Cl )-0 2 ( Cl ) = 0. (A6) 

Finally, we will consider the function 

d 

ip 3 (c) = c-— X H(c). (A7) 

From the equation obeyed by ^(Aci), we can take derivate with respect to A, and evaluate the result for A = 1. We 
arrive at an equation for ^3 (ci), 

A(ci)Mci) = (d+3)!(^-) Xh{ci). (A8) 

From equations (I A4|l . (|A6[) and (|A8[) . we can identify two eigenfunctions of A. Making use of (|A4|) and l|A8p . it 
appears that 

A(c)(~.[cxh(c)]+x*(c)) =0. (A9) 

Hence, from Eqs. (|A6|) and (|A9[) we can conclude that the null eigenvalue is (d + l)-fold degenerate with the 
eigenfunctions 

Ci(c) = ~ ■ [c X h(c)] + xh(c), 6 - -§~ c Xh{c). (A10) 
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Appendix B: EVALUATION OF THE COEFFICIENT a 2 , i2 ,i 

In this appendix, we show (£ 2 ,i(ci)£ 2) i(c 2 )|r(ci, C2)) = 0. The integral corresponding to the second term of T is 
simply 

f f ~ d d ~ 

J dc x J dc 2 c M c 2ii £ — • -^Xh{ci)xh{c 2 ) = Co- ( B1 ) 

The other term can be written as 

J dci J dc 2 c 1 .iC2.if (c ll c 2 )xH(ci)xH(c2) 

= J dc-i J dc 2 XH(ci)XH(c2) J dd-6(o- • ci 2 )(ct • C12) [6* - 1] c M c 2;i 
= J dc 1 J dc 2 XH{ci)XH{c2) J dfr<d(fr ■ c 12 ) 



(cr • c 12 ) — — (<X • c 12 )c l2ti (j l 



n (d-l)/ 2 1 _ a2 
r(^) 4rf 



1 (o- • C12) CT, 

tfci / dc 2 XH{ci)XH{c 2 )c\ 2 =g, (B2) 



which is the desired result. 



Appendix C: EVALUATION OF THE COEFFICIENTS a^, 

In this appendix we evaluate the coefficients appr. As the number of particles and the total momentum are conserved 
quantities in our system, we have 

((6N) 2 ) = 0, (SPiSPj) = 0, (CI) 

(5N5Pi) = 0, (SNSE) = 0, (C2) 

(SESPi) = 0. (C3) 

Enforcing the above constraints, we obtain 

JdcxH(c)+ J dci J dc 2 (j) H (ci,c 2 ) = l + a n = 0, (C4) 

JdcaxH(c)+ j ' dci f dc2C U <f>H(ci,c 2 ) = a 12 = 0, (C5) 

J dcc 2 XH(c) + J dci J dc 2 cj(j) H (ci,c 2 ) = ^ + dai 3 + ^an = 0, (C6) 

y dcaCjXH(c) + j ' dci J dc 2 ciiC 2j <j} H (ci,c 2 ) = ^S i3 ■+ a 2l2j = 0, (C7) 

J dcciC 2 XH(c) + J ' dci J dc 2 ciiC20ff(ci,c 2 ) = a 23 = 0. (C8) 
As a consequence, the values of some coefficients follow 

an = -1, ai2 = 0, ai 3 = --, a 2i2 j = a 23 = 0. (C9) 

Of course, the coefficients associated to A^ + A^ 7^ could also have been calculated directly by equation ([57]). 
obtaining the same results. The coefficient 033 is evaluated using (|5T[) and it can be written in terms of the one- 
particle distribution function as 

(C3(ci)6(c 2 )|To(ci,c 2 )xh(ci)xh(c 2 )) 
° 33 ~ 3Co 
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where 

d-l 



1 

with 



K a ) = ~]T p + 5N d2 J dc l J dc 2XH(ci)XH{c2)${c 1 ,C 2 ), (Cll) 



. (l-a 2 )(d+l + 2a 2 ) 5 
<?(ci,c 2 ) = ^ -4 

(d + 5)-a 2 (d+l)+4a^ 3 



-c 12 C 



2 



(2d + 3-3a)c 12 (C-c 12 )", (C12) 

and C = (ci + c 2 )/2. The coefficient 6(a) can be evaluated using the expression of xh{c) in the first Sonine approxi- 
mation, Eq. (fTT|) , which yields 

(1 + rf)(3 + rf)(2d(a 2 + 16(-1 + a) + 15a 2 a) + 16(-1 + a)(-l + 2a 2 ) + a 2 (7 + 3a(-13 + 10(-1 + a)a)))r[(l + d)/2] 
{a) ~ 2 1 /2 7r d+i/2i28d2(_ 2 + (5 + d)/2)(-l + (5 + d)/2)r[d/2]r[-2 + (5 + d)/2] { 

If we take into account the explicit form of Co and 02, given in Eqs. (|20| and (|35p respectively, we obtain after some 
algebra 

_ -15 + 7d + Ud 2 - 3(-9 + d(9 + 2d))a + 30(1 + d)a 2 - 6(9 + d)a 3 ) 
° 33 ~ 9d(-19 + 2d(-7 + 3a) + 3a(9 + 2(-l + a)a)) ' ( ' 
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